#include "../lib/lgd.h"
#include "ctime"
#include "random"
#include "iostream"

#define M 80
#define N 100

//返回范围内的随机浮点值 注意调取函数之前要调用srand(time(0));
double random_double(double l, double t)
{
	return (t-l)*rand()*1.0/RAND_MAX + l;
}

//返回范围内的随机整数 注意调取函数之前要调用srand(time(0));
int random_int(int small, int big)
{
	return (rand() % (big - small)) + small;
}

// 普通二维数组做核矩阵
double **kernel;
double *obs;

int m_count = 0;

// 计算目标函数值
lgd_float evaluate(void *instance, const lgd_float *x, lgd_float *g, const int n_size, const int m)
{
	double sum = 0.0;
	for (int j = 0; j < n_size; j++)
	{
		sum += kernel[m][j] * x[j];
	}
	sum -= obs[m];

	for (int j = 0; j < n_size; j++)
	{
		g[j] = 2.0*kernel[m][j] * sum;
	}

	return sum*sum;
}

//定义共轭梯度监控函数
int progress(void *instance, lgd_float best_fx, const lgd_float *best_x, 
	const lgd_para *param, const int n_size, const int curr_k, const int best_k)
{
	std::clog << "fight time: " << curr_k << ", fx: " << best_fx << "\r";
	//std::clog << "\033[1A\033[K";
	return 0;
}

int main(int argc, char const *argv[])
{
	kernel = new double *[M];
	for (int i = 0; i < M; i++)
	{
		kernel[i] = new double [N];
	}
	obs = new double [M];

	srand(time(0));
	// 添加一些大数
	int tmp_id, tmp_size;
	double tmp_val;
	for (int i = 0; i < M; i++)
	{
		tmp_size = random_int(25, 35);
		for (int j = 0; j < tmp_size; j++)
		{
			tmp_id = random_int(0, N);
			tmp_val = random_double(-10, 10);

			kernel[i][tmp_id] = tmp_val;
		}
	}

	// 生成一组正演解 包含一些大值和一些小值
	double *fm = new double [N];
	int N2 = (int) N/2;
	for (int i = 0; i < N2; i++)
	{
		fm[i] = random_double(5, 10);
	}

	for (int i = N2; i < N; i++)
	{
		fm[i] = random_double(1, 2);
	}

	// 计算正演值
	for (int i = 0; i < M; i++)
	{
		obs[i] = 0.0;
		for (int j = 0; j < N; j++)
		{
			obs[i] += kernel[i][j]*fm[j];
		}
		// 添加噪声
		obs[i] += random_double(-1e-3, 1e-3);
	}

	/********************准备工作完成************************/
	// 声明一组解
	double *m = new double [N];
	double *mean_m = new double [N];
	double *stddev_m = new double [N];
	double *low = new double [N];
	double *hig = new double [N];
	for (int i = 0; i < N; i++)
	{
		m[i] = 0.0;
		mean_m[i] = 0.0;
		stddev_m[i] = 0.0;
		low[i] = 1.0; // 对解的范围进行约束
		hig[i] = 100.0;
	}

	for (int i = 0; i < N2; i++)
	{
		low[i] = 5.0; // 对解的范围进行约束
		hig[i] = 10.0;
	}

	for (int i = N2; i < N; i++)
	{
		low[i] = 1.0;
		hig[i] = 2.0;
	}

	lgd_para my_para = lgd_default_parameters();
	my_para.flight_times = 20000;

	lgd_float fx;
	int ret = slgd_solver(NULL, mean_m, stddev_m, evaluate, progress, &fx, m, N, M, &my_para, NULL, low, hig);
	std::clog << std::endl; lgd_error_str(ret);

	for (int i = 0; i < N; i++)
	{
		std::cout << fm[i] << " " << mean_m[i] << " " << stddev_m[i] << std::endl;
	}

	for (int i = 0; i < M; i++)
	{
		delete[] kernel[i];
	}
	delete[] kernel;
	delete[] obs;
	delete[] fm;
	delete[] m;
	delete[] mean_m;
	delete[] stddev_m;
	delete[] low;
	delete[] hig;
	return 0;
}